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Abstract 

This paper presents a detailed discussion of Markov parameters in system identification. 
Different forms of input-output representation of linear discrete-time systems are reviewed and 
discussed. Interpretation of sampled response data as Markov parameters is presented. Relations 
between the state-space model and particular linear difference models via the Markov parameters 
are formulated. A generalization of Markov parameters to observer and Kalman filter Markov 
parameters for system representation and identification is explained. These extended Markov 
parameters play an important role in providing not only a state-space realization, but also an 
observer/Kalman filter for the system of interest 


1. Introduction 


The basic objective of system identification is to develop a mathematical model of a physical 
system by observing its input-output relationship. A mathematical model of a system is required 
for many purposes, e.g., system analysis or controller design. Often in practice, one concerns 
with linear models since many physical systems can be described by linear or approximately linear 
equations. For a given linear system, there are several models that can be used to describe the 
input-output characteristics of that system. Depending on intended application, one type of model 
may be preferred over another. Since different types of models all describe the same system, they 
must be related in a fundamental way. The relationship between different model structures can be 
described in terms of the Markov parameters, which can be used to describe the system input- 
output map in various ways. Despite the mathematical connotation attached to its name, the 
Markov parameters are simply the pulse response functions of a discrete system to a unit pulse 
input. This paper is an overview of Markov parameters commonly encountered in system 
identification theory. 

The outline of this paper is as follows. The paper starts with a review of different ways to 
describe the input-output relationship of a linear system. This is standard material, which can be 
found in various texts . 1 - 4 The relationship between the Markov parameters and these different 
descriptions is discussed. The observer equation can also be used to describe the system input- 
output relation. The relationship between the Markov parameters of an observer and the Markov 
parameters of the system is presented. This connection is particularly important since it serves as 
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the basic foundation for several recently developed system identification methods, 511 including the 
Observer/Kalman filter Identification algorithm (OKID). The Markov parameters that are 
discussed in this paper are the discrete-time Markov parameters, but physical dynamic systems are 
described by continuous-time differential equations. The basic relationship between discrete and 
continuous-time models is described. The process of converting the original second-order 
continuous-time dynamic equation to a first-order continuous-time state space equation, and then to 
a first-order discrete-time state space equation is reviewed. The reverse process of converting a 
discrete-time identified model to a continuous-time state space representation is also described. 
The Toeplitz and Hankel matrices which are the two important matrices of Markov parameters are 
explained. The problem of realization from the Markov parameters is then described with the 
Eigensystem Realization Algorithm (ERA). 1215 Many realization techniques require knowledge of 
the Markov parameters as a starting point. Different ways of determining the Markov parameters 
or parameters that have the same structure as those of Markov parameters are presented. 518 This 
process requires either a simple interpretation of the sampled response under different special test 
inputs, or a more substantial computation procedure to obtain the Markov parameters from 
experimental data. In connection with this development, particular stochastic linear difference 
equations for identification of state space models are also presented and discussed. Basic relations 
between various stochastic models and the Markov parameters describing a linear system and the 
associated observer/Kalman filter are formulated 


2. Input-Output Representation 


Dynamic systems are modeled by continuous-time or discrete-time equations. In continuous- 
time models, the input-output relations are described by differential equations. In discrete-time 
models, they are described by difference equations. A close approximation of a continuous-time 
model can be obtained by a discrete one provided that the sampling rate is sufficiently high. A 
linear discrete system is commonly described by an n-th order difference equation, the weighting 
sequence, or a discrete state space model. Relations between the Markov parameters of a system 
and these representations are described. 

2.1 Linear Difference Equation 

A /i-th order difference equation representation of a linear discrete time-invariant system with q 
outputs and m inputs takes the general form 

y(k) + Aiy(k-l)+ ••• + A m y(k-n) = B 0 u(k) + B l u(k-l) + - +B„u(lc-n) 


or 

y(k) + ^ A,y(k - i) = £ B,u(k - i ) 


where the input variable is denoted by u(k), the output variable by y(k), k is the time index, and 
Ai, B, are constant matrix coefficients of dimensions qxq, and qxm, respectively. If a delay 
operator operating on a variable z(k) is defined as 

q~'z(k) = z(k- 1) 

then Eq. (1) can be written in the following form 
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k(q~ x )y(k) = B(q l )u(Jc) 


( 2 ) 


where the polynomials of the delay operators are given as 

A(< 7 ~') = / + A x q~' + A 2 q~ 2 + + A.?”" 

B {q ’) = Bo + B^q 1 + B%q 2 + ••• + B„q 

This model of a linear discrete system is sometimes referred to as a Deterministic Auto-Regressive 
Moving Average model (DARMA). 4 This form is commonly used in developing recursive system 
identification techniques. The above description bears close resemblance to the transfer function 
representation of a linear discrete system. The transfer function H(z) is defined to be such that 
Y(z) = H(z)U(z), where U(z) and Y(z) denote the Z-transform of the input and output sequences, 
respectively. FromEq. (1), 

H(z) = (/ + Az~' + - +A n z- n y'(Bo + B 1 z- 1 + .» +B h z~ h ) (3) 

In the single-variable case, A, B, are scalar elements, denoted by a,, b,, the transfer function H(z) 
is simply 


H( _, Y(z) r , ^+^ 1 z' 1 + 

U(z) 1 + fljZ -1 + ••■ +fl,z _ * 


2.2 The Weighting Sequence Description 

For continuous-time systems governed by linear differential equations, if the response of the 
system to a Dirac delta impulse is known, its response to any other general input can be determined 
by the convolution integral. Similarly, for discrete-time systems, analogous to the Dirac impulse 
function is the Kronecker unit pulse sequence defined as 


S(k) = 



, k = 0 
, k ^ 0 


(4) 


If the system responses to a pulse sequence applied at each input are assembled to form the unit 
pulse response sequence denoted by 


K(0), T(l), 7(2), ... 

then the input-output relationship can be expressed in terms of this pulse response sequence. 
Assuming zero initial conditions, this is given by 

y(k) = Y (0)u(k) + Y (l)u(k - 1) + 7 (2)u(k - 2) + ... 

= jry(i>(*-0 (5) 
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It can be seen from Eq. (5) that the contribution to the output at time step k by the input applied at 
time step k and at previous time steps k — 1, k — 2, ... are weighted by the pulse response 
sequence. For this reason the pulse response sequence is also known as the weighting sequence, 
and this input-output description is called the weighting sequence description. If the system is 
asymptotically stable, then the infinite summation in Eq. (5) can be approximated by a finite one 
since the weighting sequence may be truncated after a finite number of time steps. For lightly 
damped systems, however, the retained portion of the weighting sequence can be excessively 
large, and this description becomes cumbersome. Taking the Z-transform of both sides of Eq. (5) 
assuming zero initial conditions yields 

Y(z) = [r(0) + 7(l)z- 1 +r(2)z- 2 + ... lU(z) 

( 6 ) 

= H(z)U(z) 


Thus the Z-transform of the weighting sequence is precisely the transfer function of the system in 
the z-domain. Furthermore, it is possible to relate the weighting sequence of a system to the 
coefficients of its difference representation in the time domain. Using Eq. (6) and Eq. (3) yields 

(/ + 4z -, + +A n z- H )(Y(0) + Y(l)z- i +Y(2)z~ 2 + ...) = B 0 + B l z~' + — +B,z _ " 


By comparing like powers of z 1 , the following relationship between the weighting sequence and 
the coefficients of the difference equation representation in Eq. (1) is obtained. 




b, , 
o , 


t = 0, 1, 2, ..., n 
i > n 


(7) 


A 0 = I 


The quantities Y(k), k = 0, 1, 2, ... are called the Markov parameters of the system. It 
immediately follows from the above discussion that (1) the Markov parameters are precisely the 
values of the weighting sequence at different time steps, (2) the sequence of Markov parameters is 
the pulse response of the system to a unit pulse, (3) the Z-transform of the Markov parameter 
sequence is the system transfer function in the complex z-domain. 

2.3 State-Space Representation 

The input-output relationship can also be expressed 
order difference equations of the form 

jc(* + 1) = Ac(*) + 

y{k) = Cx(k) + 

where the dimensions of A, B, C, and Darenxn, nxm, 
for the output y(k) in terms of the previous inputs yields 

k 

y(k) = ^CA‘- l Bu(k-i) 

i=l 


in terms of a set of n simultaneous first 


Bu(k) 

Du(k) 


( 8 ) 


qxn, and qxm, respectively. Solving 


+ Dw(Jfc) 


(9) 
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The Markov parameter sequence is simply the pulse response of the system, therefore, they must 
be unique for a given system. Comparing Eq. (9) to Eq. (5) and noting that u(k) = 0 for k < 0 
yields 


y(0) = D 

Y(k) = CA k ' l B , k = 1,2,3,... 


( 10 ) 


Equations (10) express the Markov parameters in terms of the system discrete state space matrices 
(Ay B, C, D). The state space matrices, however, are not unique since the state vector is 
coordinate-dependent. For example, let the state vector be transformed by a coordinate 
transformation T 


z(k) = Tx(k) 


Then Eq. (8) can be written as 


z(k + 1) = TAT~ l z(k) + TBu(k ) 
y(k) = CT~ l z(k) + Du(k) 


which is a state space description (TAT~\ TB, CT~\ £>) relating u(k) and y(k) by a new state 
vector z(k). The system Markov parameters computed using these state space equations are the 
same as before, 

Y(k) = CT" 1 (TAT' 1 f~ l TB = CA k ~ l B 


Thus, there are an infinite number of state-space representations that produce the same input-output 
mapping given in Eq. (9). Of particular interest are the canonical forms where the Markov 
parameters appear explicitly in these representations. For example, consider the single variable 
case, where the transfer function H(z) in Eq. (3) is rewritten as 


z" + a,z" + + oc H -\Z + oc H 

The observability canonical form of the above system is 



i 

o 

o 

« 

• 

• 

o 
1 


T(l)' 


0 0 1—0 


Y( 2) 

A) — 

0 0 0 : 
• * • • M 

I • • • B x 

i Bo — 

n 3) 


-a H -a,- 2 — ~oci_ 


/(»)_ 


Co = [l 0 0 — 0] Do — d 


( 11 ) 


( 12 ) 


where the coefficients in the input matrix are precisely the first n Markov parameters obtained from 
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( 13 ) 


gifll+i +P* z iz + P 1 _ y r(j . )2 
z" + a 1 z '" 1 + +a„_,z + a:„ " 

However, in the multiple variable case, such a canonical representation with the Markov 
parameters appearing explicitly is not of minimum order. 

It is also possible to see the Markov parameters in the transfer function description of the 
system in state space format as given in Eq. (8). Taking the Z-transform of Eq. (8) yields 
Y(z) = H(z)U(z) where the transfer function H(z) can be expressed as 

H(z) = D + C(zI-AY l B 

. ' 2 3 (14) 

= D + CBz~ l + CABz " 2 + CA 2 Bz ~ 3 + - 

Finally, consider the case where the system is asymptotically stable. In this case, 


limC4*fl = 0 

The input-output description in Eq. (9) can be approximated by a finite number of Markov 
parameters 


y(k) - Du{k) + CBu(k - 1) + CABu(k - 2) + ••• +CA P x Bu(k — p) 

= ^ CA‘~ l Bu(k - i) + Du(k) (15) 


where p is sufficiently large such that CA k B~0, k>p. By comparing with Eq. (1), this 
description is a special case of the linear difference description with 

A t =0 , B q = D 

(16) 

B i = Y(i) = CA i ~ 1 B , i = l, 2, ..., p 


The Bi coefficients appear explicitly as the system Markov parameters. Note that this description 

with a finite number of system Markov parameters is possible only for asymptotically stable 
systems. 

2.4 Observer Equation for Input-Output Representation 

The system input-output relation can also be described by a set of state space equations derived 
as follows. Add and subtract the term My(i) to the right hand side of the state equation in Eq. (8) 
to yield 

x(k + 1) = Ax(&) + Bu{k) + My(k) - My(k) 

= (A + MC)x{k) + (B + MD)u(k) - My(k) (17) 

Define the following matrices 
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(18) 


A=A+MC 


B=[B + MD, -M] 



Then the original system becomes 

x(k + 1) = Ax(k) + Bv(k) 

y(k) = Cx(k) + Du(k) 


(19) 


If M is a matrix such that A = A + MC is asymptotically stable then the input-output relation can be 
approximated by a finite number of parameters CA ‘~ l B , say p of them. 

y(k) = £ CA‘- ] Bv(k - i) + Du(k) (20) 

»=1 


Noting that v(k) contains both u(k) and y(k), Eq. (20) can be rewritten as 

y(k) + ^C(A + MCj~ l My(k - i) = ^C(A + MCj~\B + MD)u(k - 1) + Du(k ) (21) 

;=i >= l 


Observe that Eq. (20) is in a linear difference form as in Eq. (1) with the coefficients given by 


Ai = C(A + A/C) 1 " 1 M , B 0 = D 
Bi - C(A + MC)" 1 ( B + MD ) , / = 1, 2, .... p 


( 22 ) 


The role of the matrix M in the above development can be interpreted in terms of an observer. 
Consider the system given in Eq. (8). It has an observer of the form 


x(k + 1) = Ax(k) + Bu(k) - M[y(k) - y(*)] 
y(k) = Cx(k) + Du(k) 


(23) 


The state equation in Eq. (23) can be expressed as 

x(k + 1) = Ax(Ar) + Bu(k) - MC[x(k) - x(*)] 

= (A + MC)x(k) + Bu(k) - M[y(k) - Du(k )] (24) 

= (A + MC)x(k) + (B + MD)u(k) - My(k) 


Defining the state estimation error e(k) = x(k)-x(k), the equation that governs e(k) is 

e(k + 1) = Ax(k) + Bu(k) - [(/4 + MC)x(k) + (B + MD)u{k) - My(k )] 

= (A + MC)e(k) 
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If the system (8) is observable, then the matrix M may be chosen to place the eigenvalues of 
A + MC in any desired (symmetric) configuration. If the matrix M is such that A + MC i s 
asymptotically stable, then the estimated state x(k) tends to the true state x(k) as k tends to infinity. 
Equation (24) then becomes 

x(k + 1) = ( A + MC)x(k) + (B + MD)u(k) - My(k) 

which is exactly the same as Eq. (17). From this analysis, the matrix M can be interpreted as an 
observer gain. 

The parameters CA'~ X B in Eq. (19) are the Markov parameters of an observer system, hence 
they are referred to as observer Markov parameters. They are denoted by 

7(0) = D , Y(k) = CA k -'B , k = \, 2, p (25) 

The extra freedom inherent in these parameters can be exploited to develop various identification 
algorithms. 5-11 Consider the case where M corresponds to a deadbeat observer gain, the observer 
Markov parameters become identically zero after a finite number of terms. For lightly damped 
structures, this means that the system can be described by a reduced number of observer Markov 
parameters Y ( k ) instead of an otherwise large number of the usual system Markov parameters 
Y(k). For this reason, the observer Markov parameters are useful in the development of various 
system identification algorithms. The relationship between the Maikov parameters of the observer 
system and those of the actual system will be presented in the following section. 


3. Relationship between the Markov Parameters 
of the Observer and Actual System 


The observer Markov parameters contain information about the system and the observer. It is 
possible to recover both the system and the observer gain from a given set of observer Markov 
parameters. This section describes this relationship. 

3.1 Recovery of System Markov Parameters from Observer Markov Parameters 

In this section, the basic relationship between the Markov parameters of the observer and 
actual system is described. First, it is shown in the following that given a set of observer Markov 
parameters, the system Markov parameters can be recovered uniquely. The direct transmission 
term D is simply 


D = Y(0) = Y(0) 


(26) 


For ease of presentation, the observer Markov parameters are partitioned as 


Y(k) = CA k -’B 

= [C(A + MC) k ~ 1 (B + MD) -C(A + MC) k ~ l M\ 

= [F (1) 0t) 7 (2) (k)] (27) 
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The Markov parameter 7(1) = CB of the system is simply 

y(l) = CB = F (1) (l) + F (2) (1)7(0) 

To obtain the next Markov parameter 7( 2) = CAB, consider Y (2) 

7(2) = CAB = 7 (1) (2) + F (2> (1)7 (1) + F <2) (2)T (0) 


Similarly, 


7( 3) = CA 2 B 

= F (1) (3) - CMCAB - C(A + MC)MCB - C(A + MCf MD 
= F (1) (3) + F (2) (1)7(2) + F (2) (2)7(1) + F (2) (3)7(0) 

By induction, the general relationship between the actual system Markov parameters and the 
observer Markov parameters can be shown to be 

Y(k) = Y m (k) + 2 7 <2) (i)Y(k - i ) (28) 

;=i 

where Y(k) are considered to be zero for k > p. The above recursive equation can be written in 
matrix form as 


" I 

"7(1)' 


'F (1) ( 1)' 


'F (2) (l)‘ 

-F (2) (l) / 

7(2) 


F (1> (2) 


F (2) (2) 

— F (2) (2) -F (2) (l) 

* • ' r 

7(3) 

= 

F (,) (3) 


F (2) (3) 

1 I i 

-F (2) (jfc-1) -Y m (k- 2) ••• -F <2) (1) I 

Y(k\ 


_F (1) (*)_ 


, 

w 

1 


The left most matrix in the above equation is square and full rank, which implies that from a given 
set of observer Markov parameters 7 (k), the system Markov parameters 7(k) can be uniquely 
recovered. 


Furthermore, let the matrices Y (2) , and Y be defined as follows 

Y (2) =[-F (2) (p) -F (2) (p-1) -F (2) (p-2) — -F (2) (l)] 


Y = [7(p + 2) 7(p + 3) 7(p + 4) Y(2p + \)] 


then from Eq. (28) 


Y (2) H(1) = Y 


(30) 


(31) 


where 
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' Y( 2) 

7(3) 

7(p + l) 

H(l) = 

7(3) 

7(4) 

m 

w 

- 7(p + 2) 


Y(p + 1) 

7(p + 2) 

- 7(2p) 


( 32 ) 


Examining the rank properties of Eq. (31) reveals that to describe a system of order n, the number 
of observer Markov parameters p must be such that qp^n where q is the number of outputs. 
Furthermore, the maximum order of a system that can be described with p observer Markov 
parameters is qp. The implication of this analysis is that for multiple-output systems, the number 
of observer Markov parameters required to be identified can be substantially less than the true order 
of the system. Specifically, the minimum number of observer Maikov parameters that can describe 
the system is which is the smallest value of p such that qp^ > n. 

3.2 Recovery of Observer Gain from Observer Markov Parameters 

Since the observer Markov parameters contain information of both the system and observer 
gain, it is possible also to recover the observer gain M from the observer Markov parameters. 
First, the sequence 


Y M (k) = CA k ~ l M , * = 1,2,... 


(33) 


can be computed from the observer Markov parameter sequence Y (k) = CA k l B, k = 1, 2, ... as 
follows 


Y„(k) = ~Y m (k) + i £Y m (i)Y„(k - i) (34) 


The parameters in Eq. (33) contain information about the observer gain, therefore they are referred 
to as observer gain Markov parameters. Given A, C, and a sufficient number of observer gain 
Markov parameters Y u (k) computed from Eq. (34), the observer gain M can be obtained from 

M = (V k T V k y'v k T Y M (35) 


where 



" C ' 


' CM ' 


CA 


CAM 

v* = 

CA k 

Y„ = 

CA k M 


provided V k is of rank n. Alternatively, the observer gain matrix M can be realized simultaneously 
with the system matrices A, B,C from the combined Markov parameters 
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( 36 ) 


P(k) = [CA k -'B CA k ~ l M] 

= CA k -'[B M] , k = l, 2, ... 


This problem will be revisited in a later section. 


4. Relationship between Discrete and Continuous-Time Models 


The results presented so far are concerned with discrete-time models. However, all structural 
systems are continuous, but information available for identification is in the form of sampled data. 
In this section, the relationship between continuous and discrete-time models is described. 

4.1 Continuous-Time Dynamic Equations 

Consider the equations of motion for a finite-dimensional linear dynamic system of the form 

Mw(r) + Dvv(r) + Kw(/) = Bfit(t) (37) 


where M, D, and K denote the system mass, damping, and stiffness matrices, respectively. If the 
dynamic system is measured by output quantities such as position, velocity, and acceleration then 
in general the set of observation variables can be written in matrix form as 

y(0 = CMt ) + C 2 w(t) + C 3 w(t) 

= [C, - C 3 M -1 K]w(0 + [C 2 - C 3 M _l D]w(f) + C 3 MT'B f u(t) 

Equation (37) together with the ouput equation can be written in first order form as 


where 


x(t) = A c x(t) + B c it{t) 
y(t) = Cx(t) + Duit) 


x(t) = 


w(r)' 

vHoJ ’ 


■ 0 

/ ' 

. Br = 

" 0 ' 

-M -1 K 

-M"'D 


M~'B f 


C = [C,-C 3 M-'K C 2 -C 3 M-'D] , D = CM~'B f 


(38) 


(39) 


This is the continuous-time equations in state space format. The continuous-time state space 
equation can be converted to the discrete-time equation by the following procedure. 

4.2 Conversion of Continuous to Discrete-Time State Space Model 

The general solution for x(t) at any time t for any input u(t) from any initial condition x(t 0 ) is 


1 1 


X(0 = e A ‘ (, " a) x(to) + f e^'-^BM^dr 
J «8 

Consider discrete sampling time intervals 0, 7, 27, .... To see how x(t) changes from one time 
step to the next, replace t 0 = kT and t = (k + 1 )T in the above equation. Furthermore, assume a 
zero-order hold on the input, i.e., the input u(kT) is held constant over the time interval from kT to 
(k + 1)7 , one obtains 

x[(k + 1)7] = e^xikT) + e Mii+l)T ~ x] B c u(t)dT 
= e KT x(kT) + ^e^dT'B c ) ju(kT) 

Using a simplified notation k for the time argument kT , 

x(k + l) = Ax(k) + Bu(k) 
y(k) = Cx(k) + Du(k) 

where the output equation is written at every time step, and 

A = e KT , B = \ T e^dxB, 

JO 


(40) 


(41) 


and the C and D matrices are the same as in the continuous-time representation, Eq. (38). Equation 
(40) is a discrete-time state space equation with a sampling interval 7 for the continuous-time 
dynamic system in Eq. (37). 

In practice, to compute the matrix exponential, one can use an approximation of the form 19 


e 


KT _ 



y>o 


(42) 


where Rm(z) is the (d.d) Padi approximant to e l 





(2d-fc)ld! 

(2t?)!fc!(t?-&)! 


(43) 


In fact, the discrete-time A and B can be computed from the continuous-time 4: and B c in one step 
by making use of the following identity 




<D (r = 7) = e L ° 



(44) 
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The above identity can be proved by first noting that <t>(0 is the solution of the following matrix 
differential equation 


JO 

dt 


A c B c 
0 0 


o 


0(0) = / 


(45) 


Denote O n , 0 12 , 0 2I , O^ as the partitions of O, 

O = 


o„ o 12 

0 21 o 


From Eq. (44) 


JOji . _ jo 12 . . _ . jo 21 

— — = A c O u , —^ = 40,2 + 5,022 , — — 

dt dt dt 


= 0, 


JO 22 _ 


dt 


= 0 


(46) 


Application of the initial conditions O u = /, 0 21 = 0, O^ = I yields On = A, 0 2 , = 0, 0 22 = I . 
To show 0 12 = B, pre-multiply both sides of the 0, 2 equation by e~ Ac> and rearrange terms 

-A, JO,2 -A t t A ^ __ -a c , 


dt 


-^‘40,2=^ 


Integrating both sides of the above equation from t = 0 to t = T, recognizing that the left hand side 
is a perfect differential, yields 


m w]<* - = iy jB ‘ dt 


<=T rT 


dt 


Solve for 0,2 (T), noting that O 12 (0) = 0, 

0 12 (T) = j \ a ' (t - ,] B c dt = jV‘ T 5 c dt=B 


(47) 


The last step in Eq. (47) is obtained by a change of integration variable from (T - 1) to T. 

4.3 Conversion of Discrete to Continuous-Time State Space Model 

A realized discrete-time system matrix A can be converted to a continuous-time one. Assume 
that A can be diagonalized by a matrix P such that P~ l AP = A = diag(X u A 2 , A*), then 

A = P~ l e A < T P = e PA ' PT =e [t t t} 

A continuous-time system matrix 4 can be obtained from 
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Ac = Pdiag 


InAi 

7 ’ 


lnA 2 
7 ’ 



The eigenvalues of the continuous-time A arc 

, . In A, 
y J = o J ± i(o j = —— 




2 In 


1 = 0 , 1 , 2 , 


(48) 


(49) 


There is no unique continuous-time system that produces a given discrete-time model. Under the 
assumption that the sampling rate is sufficiently fast such that all frequencies in the signals being 
identified fall below the Nyquist frequency, an appropriate value of l may be chosen. In practice, 
there will always be frequencies outside the Nyquist frequency, and this choice of l will 
misinterpret as frequencies within the Nyquist frequency and thus produce modelling errors. This 
phenomenon is known as aliasing. In the identification of structures, the complex parts of the 
eigenvalues are the modal frequencies, and the real parts are the modal damping factors 
characterizing the decay rates of the modal oscillations. 


To convert the discrete B to a continuous B c , first write A = e A '‘ as a convergent power series 


A = e A ‘‘ = I + Aj + 


Ay | Alt 3 
2! + 3! 


(50) 


Therefore, 


Hence 


4 . , „ A 2 e T 2 A 3 c T 3 

A- 1 = A C T + — — + — — + 
2! 3! 


a 2 b t 2 a 3 r t 3 

(A - I)B C = A e B c T + + .. + 


2! 


3! 


(51) 


Integrate the second equation in Eq. (41) term by term, and pre-multiply by A c 


A C B — A C B C T + 


AlBcT_ f AMI [ 

2 3x2! 


Comparing Eq. (51) and Eq. (52) immediately yields 

(A-I)B c = AcB 

The continuous-time input matrix B c is 


(52) 


B C = (A — iy 1 AcB 


(53) 


14 



In this section, to derive Eq. (48) and Eq. (53), the matrix A is required to be diagonalizable. This 
requirement excludes the presence of rigid body modes in the model. To account for rigid body 
modes, a different procedure must be used to convert from a discrete to a continuous-time 
presentation. Note that the reverse procedure to convert from the continuous to a discrete-time 
presentation formulated in the previous section using matrix exponential does not have this 
limitation. 


5. Toeplitz and Hankel Matrices 

Two important matrices often used in system identification and control are the Toeplitz matrix 
and the Hankel matrix. They are simply matrices of Markov parameters of particular structures. 

5.1 The Toeplitz Matrix 

Often associated with the Markov parameters is a lower block triangular Toeplitz matrix. The 
input and output histories of the discrete system given Eq. (8) satisfy the relation 


'y(k) 


' C ' 


' D 

’u(k) 

y(k + 1) 

- 

CA 

x(k) + 

CB D 

9 • • 

* • • 

u(k + 1) 

_y(k + n- 1)_ 


_CA m ~ l 


CA h ~ 2 B CB D 

_u(k + a — 1)_ 


The observability matrix V H . X and the Toeplitz matrix T(n-l) formed by the Markov parameters 
are defined as 



‘ C ' 


' T(0) 

V_i* 

CA 

T(n-1) = 

T(l) f(0) 

9 9 * 

• 9 9 


CA H ~ l 


Y(n-l) .» T( 1) Y(0) 


(55) 


with the blank spaces denoting zero elements. The Toeplitz matrix is also known as the pulse 
response matrix. 

5.2 The Hankel Matrix 


Another important matrix of Markov parameters which is often encountered in realization 
theory is the Hankel matrix. The Hankel matrix is defined as 


H0fc-1) 


'Y(k) 

Y(k + 1) 


T(/t + l) 
Y(k + 2) 


— Y(k + s) 

- y(*+s+i) 


Y(k + r) Y(k + r + 1) 


Y(k + r + s) 


(56) 


1 5 



If r + 1 > n and s + 1 > n , the Hankel matrix is of rank n. The Hankel matrix can be factorized in 
terms of the system observability matrix V r and controllability matrix W, as 

H(k) = V r A k W, (57) 


where 



W S =[B AB — A'B] 


For the scalar case considered in the previous section, it is interesting to point out a 
relationship between the system parameters, the Hankel matrix, and the Toeplitz matrix for the 
system characterized by the transfer function given in Eq. (11). They are related by 


H(0 )[a„ a„_, tt]] r = [Y(n + l) Y(n + 2) - 7(2n)f 

[A A A] r = T(n)[i Qfi ■■■ a n . i] r 


(58) 


where 



'7(1) 

7(2) 

... Y(n) ~ 


’7(1) 

H(0) = 

7(2) 

7(3) 

••• 7(n + l) 

T(n) = 

7(2) 7(1) 


Y(n) 

7(n + l) 

... 7(2n — 1)_ 


Y(n) - 7(2) 7(1) 


6. Realization from Markov Parameters 


The problem of finding a state space model from a given sequence of Markov parameters is 
known as realization. The theory of minimal realization finds such a model with minimum 
dimensions. A realization technique applicable to identification of structures is the Eigensystem 
Realization Algorithm (ERA) which is summarized in the following. 

6.1 State Space System Realization by ERA 

ERA considers the problem where a sequence of Markov parameters y(l), Y( 2), 7(3), ... is 
given. The objective is to find a state-space realization A, B, C such that the relation given in Eq. 
(10) holds. The algorithm begins with a Hankel matrix defined as in Eq. (56). The order of the 
system n is determined from the singular value decomposition of H(0) 

H(0) = U'LV T =U x Yj X Vy (59) 
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The columns of U x and V; in Eq. (59) are orthonormal, and correspond to the positive singular 
values in an n x n diagonal matrix Si . A discrete-time minimal order realization of the system can 
be shown to be 


K = Z? I2 U? mWx Sr' /2 , B r = 1} n V?E m , C r = E]U X Sj /2 (60) 

where El = [I a O a O a ], a = m, q. Thus, for a finite dimensional system, knowledge of a 
sufficient number of Markov parameters is adequate to obtain a state-space realization for the 
system. 


An important aspect in connection with system realization from Markov parameters is the order 
of the minimal realization. In the noise-free case, given a sufficient number of Markov parameters, 
the minimal order n is equal to the rank of the Hankel matrix, which is equal to the number of 
positive singular values in Si . Equation (59) holds exactly. This is evident from the factorization 
given in Eq. (57) and the fact that a minimal order realization must be observable and controllable. 
In practice, however, all structural systems have noises and nonlinearities, the problem of rank 
determination is not trivial. The singular value decomposition step in the above realization 
procedure is normally used to determine this order. Certain smaller singular values corresponding 
to less significant modes are attributed to noises, and truncated. The number of retained singular 
values then determine the system order. The corresponding retained columns of U and V in the 
singular value decomposition of H(0) are used in Eq. (60) to obtain a realization. 

6.2 State Space System/Observer Realization 

Realization can be thought of as a factorization of a parameter sequence of the form 
Y(k) = CA k ' x B, to obtain a set of (C, A, B) that preserves the prescribed relationship between the 
parameters in the sequence. Therefore, the theory is applicable to factorization of all such 
sequences whose parameters obey the same prescribed relationship. The ERA procedure can be 
applied to the combined Markov parameter sequence P(k ) given in Eq. (36) by first forming the 
Hankel matrix of P(k) 


H e (*-1) = 


P(k) 

P(/t + l) 

P(k+s) 

P(k + 1) 

P(k + 2) 

9 

9 

P(k + s + 1) 

/>(* + r) 

• 

P(k + r + 1) — 

P(k + r + s) 


(61) 


The following realization will simultaneously identify A, B, C and the observer gain M 

4 = A- 1 'W(l)aA- 1/J , [Br M r ] = Dt n QlE m+ , , C r = ElP D\ n (62) 

where H c (0) = PDQ T = P x D\Q[, and El = [/ a O a O a ], a = m + q, q. 

In the previous development the set of parameters is precisely the Markov parameter sequence 
of a discrete system. In system identification, the Markov parameters are obtained or computed 
from measurement data. The determination of the Markov parameters from sampled data is 
discussed in the following section. 
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7. Determination of Markov Parameters from Sampled Data 


Under certain conditions, it is possible to interpret the sampled data directly as Markov 
parameters. In some other cases, the Markov parameters can be computed from input-output data. 
In particular, the following cases are considered. 

7.1 Impulse Test Input 

Let the system given in Eq. (38) be subjected to an impulse in the y-th input, starting from zero 
initial conditions. The response is recorded at evenly spaced discrete-time intervals, and denoted 
by y°\kT ), k = 0, 1,2, ... beginning with y O) (0) = y U) ( 0 + ) immediately after the impulse has 
entered the system. The sampled response for the y-th input impulse is 

kT 

y u \kT) = Ce KW x(0) + Ce^ r j e~ A ‘ T B ( c Mr)dr+ Du(kT) 

0 (63) 

= C(e A ’ T ) k B ( c J) = CA k B ( J ) 


where B { C J) denotes they'-th column of the input matrix B c , j = 1, 2, ..., m. The same impulse test 
is applied to every input in the system, one at a time. The data is then assembled to form the 
following data sequence, which has the same structure as that of the Markov parameters 

[y (,) (*T) y™(kT) ••• y {m \kT)] = [CA k B? ) CA k B? ) — CA k B^ n) ] 

= CA k [B c (1) $ 2) — 

= CA k B c , * = 0, 1, 2, 3, ... (64) 

where A is the discrete-time system matrix for the system in Eq. (38) corresponding to the 
sampling period T, A = e*' 7 , but B c is the continuous- time input matrix, and C the output matrix 
which is the same for both discrete and continuous-time cases. Note that the data is recorded after 
the impulse has entered the system, therefore, information regarding D is not available. Since a 
true impulse is actually not available in practice, for a sufficiently rich input, the needed data 
sequence can be obtained by taking the inverse fast Fourier transform (FFT) of the ratio of the FFT 
transforms of the input and output measurements, 

7.2 Unit Pulse Test Input 

A way to obtain the true Markov parameter sequence is by applying a unit pulse for the y-th 
input 


Uj(kT) = 


, k = 0 

, k = 1, 2, 3, ... 


(65) 


with a zero-order hold on the input to maintain the input values during each sampling interval. The 
discrete-time equivalence of the continuous-time system of Eq. (38) for the y-th input is 



x(k + 1 ) = Ax(k) + B 0) Uj(k) 
y U) (k) = Cx(k) + D U) Uj(k) 

where B U) and D U) denote the y'-th column of B and D, respectively. The corresponding response 
is measured at discrete sampling intervals, and denoted by y U) (k.T) t k — 1, 2, 3, .... With zero 
initial conditions, a separate test is conducted for each input. The sampled responses of the 
continuous-time system are then assembled to form the data sequence 

[y (,) (0) y (2) (0) — y ( "°(0)] = [D (1) D m D {m) ] = D (66) 


and 


[y m (kT) y {2) (kT) ■■■ y (m \kT)] = [CA k ~ x B m CA k - l B m — G4* _1 fl (m) ] 

= CA l ~ 1 [B (l) B (2) ... B (m) ] 

= CA k ~ l B , k = 1, 2, 3, ... (67) 


Observe that this sequence is precisely the true Markov parameter sequence of the system discrete- 
time state space model with a sampling interval equal to the duration of the unit pulse test input. 

7.3 Extended Unit Pulse Test Input 

To put more energy into the input, consider a unit pulse in the y-th input held for p time steps 


Uj(kT) = 


, k = 0, 1, 2, .... p - 1 
, k = p, p + 1, p + 2, ... 


( 68 ) 


With zero initial conditions, the response sequence at time steps 0, T, 27, ..., (p - 1)7 is 

y O) (0) = D in 
and 

y U) ( T) = CB U) + D U) 
y U) (27) = CAB lj) + CB U) + D U) 


y U) [(p - 1)7] = CA p ~ 2 B U) + — + CB U) + D U) 


and for time steps beyond pT 


y u \pT) = CA p - l B U) + ... + CB U) 
y U) [(p + 1)7] = CA P B U) + — + CAB (i) 


The Markov parameter sequence associated with the y-th input is simply 
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D U) _ y(/)( 0) 

CB U) = y U) (T) - y U) (0) 

CA p ~ 2 B (i) = y (J) [(p - 1)7] - y U) [(p - 2)T] (69) 

and for k = p- 1, p , ... 

CA k B (i) = y U) [(k + 1)1] - ^CA k ~‘B U) (70) 

t=l 

The results for each input are then assembled to obtain the system discrete-time Markov 
parameters. 

7.4 Free Decay Response 

Consider the free response of the system given in Eq. (40) from a non-zero initial condition, 
the sampled response is 

y(hT) = CA k x( 0) , * = 0, 1, 2, ... (71) 

The above sequence has the same mathematical structure as the Markov parameter sequence, except 
that the initial condition x(0) plays the role of B. An identification procedure applied to the above 
sequence will produce a realization of C, A, and the initial condition x(0). 

7.5 Random Test Input 

Consider the case where the system given in Eq. (40) is driven by an independent, white, 
zero-mean, stationary random input sequence with unit covariance 

E{u(kT)u(kT) T } = I , k = 1, 2, 3, ... (72) 

With a zero-order hold on the input, the resultant response is measured at discrete sampling 
intervals, y(kT), k = 1,2, 3, .... First, note that 

£{y(*)«(*) T } = E{[Cx(k) + Du(k)]u(kf} = D (73) 

Next, consider 

E[y{k + 1 )u(k) T ] = £{[C4jc(*) + CBu(k) + Du(k + l)] W (*) r } = CB 

In general, the Markov parameters of the system are related to the covariance between the input and 
output data by 

£{y0t + i>(O r } = C4*- 1 5 , k = 1,2,... (74) 
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In practice, if the above relations are to be used to compute the Markov parameters, the expected 
values are replaced by sample covariance computation. This approximation is only valid when the 
input-output data sequences are stationary and the records are sufficiently long such that the 
ergodicity property of stationary random processes apply. 

7.6 System Markov Parameters from General Input-Output Data 

The problem of determining the Markov parameters from a general test input is can be 
addressed using various input-output models that are described in previous sections. Suppose that 
a set of N + 1 measurements of y(i) and u(i ), i = 0,1, N, is available. Assuming zero initial 
conditions, the input-output data can be written in matrix form as 

y = Y*_ 1 U A , (75) 

where 


y=[y(0) y(i) 

... yip) 

yip+l) ." 

y(N)] 

Y*_,=[D CB - 

.. CA p ~ l B 

CA P B - 

CA n -'b] 


’m(O) m(1) ••• 

u(p) 

u(p + 1) 

... u(N) 


u( 0) 

u(p- 1) 

u(p) 

... u(N- 1) 

u* = 


«(0) 

u( 1) 

: : 




u( 0) 

*• • 


«( 1 ) 
u( 0 ) 


In the single-input case, U* is a square matrix. The Markov parameters in can be solved 
from this equation, 


Ya,_! = yU* < 76 > 

provided w(0) *0so that U w is full rank. This solution cannot be extended to the multiple-input 
problem since in this case Y S -i cannot be determined uniquely. If the system is asymptotically 
stable, and p is sufficiently large such that the Markov parameters CA k ~ l B can be neglected for 
k> p then it is possible to solve for the Markov parameters D, CB, CAB, ..., CA P ~ X B from the 
linear set of equations 


y = YU 


(77) 


where 


Y = [£> CB — CA P ~ X B ] 
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m(0) w(1) ••• u(p ) u(p + 1) ••• u(N) 

w(0) ••• u(p-l) u(p-2) ••• M(iV-l) 

« * » * m \ 

• * • * * I 

« » v * • 

u(0) w(l) u(N-p)_ 

The solution for the system Markov parameters from Eq. (77) is 

Y = yU T (UU T ) 1 (78) 

provided that the data record is sufficiently long, and the input is sufficiently rich such that the 
inverse in Eq. (78) exists. For lightly damped, large flexible structures with multiple sensors and 
actuators, this direct solution for the Markov parameters is not desirable in practice since the 
number system Markov parameters required in Y may be excessively large. 

If the initial conditions are not zero then a slightly different equation must be used in place of 
Eq. (78) to solve for the system Markov parameters. This is given as 

Y = y,ur(u,ur ) 1 (79) 


where y, and U, are obtained by deleting the first p columns in y and U, repectively. 

7,7 Observer Markov Parameters from General Input-Output Data (OKID) 

The above mentioned problem can be overcome by solving for the observer Markov 
parameters instead. The observer equation counterpart of Eq. (77) is simply 

y = YV (80) 

where 

Y = [D CB CA pl B] 

u( 0) u(l) - u(p) u(p + Y) u(N) 

v(0) v(p-l) v(p-2) ••• v(N- 1) 

v(0) v(l) — v(N-p ) 

where v(k) is defined as in Eq. (18). A solution for the observer Markov parameters is 

Y = yV T (VV T ) -1 (81) 

If the initial conditions are not zero then a slightly different equation must be used to solve for the 
observer Markov parameters. 
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(82) 


Y = yv, r (v j vr)"' 


where y, and V, are obtained by deleting the first p columns in y and V, repectively. 

The system Markov parameters and the observer gain Markov parameters are recovered from 
the observer Markov parameters using Eq. (28) and Eq. (34), respectively. As mentioned in 
Section 6.2, application of ERA will produce a realization of both the system state space matrices 
and the corresponding observer gain. 

In the noise-free case, it is sufficient to solve for a minimum number of observer Markov 
parameters to recover the system Markov parameters. As discussed previously, for a system of 
order n, this is the minimum value Pmm such that qp min £ n, where q is the number of outputs. 
However, if more than the minimum number of observer Markov parameters are to be solved for, 
P > Pmm, then the matrix V (or V,) is row rank deficient. This implies that the observer Markov 
parameters are not unique. A solution can still be found by simply replacing the inverse of the 
quantity in the parentheses in Eq. (81) by a pseudo-inverse. 

With noises, however, the matrix V (or V,) tends to be full rank. In this case, the solution 
given in Eq. (81) or Eq. (82) is the least-squares solution that minimizes the residual 

e = y-YV or e, = y, -YV, 

depending on either Eq. (81) or Eq. (82) is used to solve for Y. In the presence of process and 
measurement noises that are white, zero- mean, and uncorrelated with the data, the least squares 
result can be interpreted in terms of the Kalman filter. The relationship between the deterministic 
and stochastic approaches is established in Ref. 10. Similar insights can be obtained by examining 
the structures of die stochastic models considered in Ref. 1 1, and in Section 8 of this paper. 

7.8 Observer Markov Parameters from General Input-Output Data with 
Prescribed Observer Poles 


The problem of identification of observer Markov parameters with prescribed observer poles is 
briefly presented here. For an observable system, it is always possible to assign the observer 
poles in any (symmetric) configuration. For simplicity, consider the case where all prescribed 
poles are real, distinct, with magnitudes less than one. They are denoted by A,, i = 1, 2, .... n. 
It can be shown that Eq. (21) can be written as 


yW = a2AL T W-/-l) + p20^-'- 1 ) + Dw W = [ a P D ] 


1=0 


T=0 


>(*- 1 )' 
<p(k- 1) 
u(k) 


(83) 


The quantities <P(k- 1) and (p(k - 1) are the weighted input and output data vectors 

<p(k - 1) = 3 m «(£ - p) , (p{k-\) = Z„y{,k-p) (84) 


where 3 m , 3, are the weighting matrices formed by the prescribed eigenvalues, 
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3/ = 


’ 2 O’- 1 ) 2 ip-2) 2 (I) r 

zLi ,t " # * 4u t / **x/ 

2 o-i) 2 (P-2) 3 (l) J 

6l 2J OLl.t "’ lk2,* *4 


/x* 


US 1 ’ is 


^ (i) / 


/X/ 


t = m, q 


X™ = diag(x k i , A?, ■•• ,A-) <x/ , i = l, 2, ..., n , k = \, 2, ..., p- 1 


and m(j — p) and y(i - p) denote the mp x land qp x 1 input and output history vectors defined as 



1 

• 1 

_ 1 


~y(k-p) 

1ST 

*- 

1 

h 

u(k- 2) 

, y(k-p) = 

y(k- 2) 


«(*-!)_ 


,y{k- 1). 


For convenience, Eq. (83) can be simply written as 

y(*) = yr(*-l) (85) 


where 


T = [o P D ] 


r(*-i) 


M- 1) 

<p(k - 1) 
u(k) 


( 86 ) 


For a set of N + 1 measurements of y(i) and u(i), i = 0, 1, .... N, assuming zero initial conditions, 
the input-output relation can be written as 


y=rr 

where y is given as in Eq. (75) and Pis 

y = [y(0) y(l) y(p) y(p + l) — y(A0] 
r=[r(-i) r( 0 ) ... r(p-i) r(p) •- w - d ] 


(87) 


(88) 


with u(k), y(k) set to zero for k < 0. The least squares solution for the observer parameters Y are 
computed from 

T=yr’-(rr r )‘ (89) 

As before, if the initial conditions are not zero then the observer parameters can still be solved from 
Eq. (89) with y and T replaced by y ( and F,, which are obtained by deleting the first p columns 
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of y and r, respectively. The observer Markov parameters are computed from the observer 
parameters as 

Y(k) = C(A + MD -M] 

= [aX<,*-'> (90) 

where 

x<r> =[Agr" iu" - AS 0 ] . <=».* 

As before, the system Markov parameters and the observer gain Markov parameters are recovered 
from the observer Markov parameters using Eq. (28) and Eq. (34), respectively. The general case 
where the prescribed poles are complex or deadbeat can be treated similarly. The readers are 
referred to Refs. 6 and 8 for further details. 


8. Stochastic Linear Difference Equations 


This section examines some particular stochastic linear difference equations for the 
identification of state space models in the presence of noises. In particular, two models are 
considered that can be interpreted in terms of a Kalman filter. 911 

8.1 Kalman Filter Model 

Consider the case where the state-space equations given in Eq. (8) are extended to include 
process and measurement noises 


x(k + l) = Ax(k) + Bu(k) + w 1 (k) 
y(k) = Cx(k) + w 2 (k ) 

The process noise Wi(k) and measurement noise w 2 (k) are two statistically independent, zero- 
mean, stationary white noise processes with covariances Q and R respectively. The same system 
can also be expressed in the form of a Kalman filter 

x (k + 1) = Ax(k) + Bu(k) + Ke(k ) 

y(k) = Cx(k) + Du(k) + e(k) (92) 


where e(k) is a white sequence of residual with covariance E = CPC 1 + R, and P is the unique 
positive definite symmetric solution of the algebraic Riccati equation 

p = apa t - apc t (cpc t + r)~'cpa t + Q 

The Kalman filter gain K is given by K = APC T (CPC T + /?)'= APC t 1T 1 . The system in Eq. 
(91) can also be expressed as 
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( 93 ) 


x(k + 1) = ( A - KC)x(k) + (B- KD)u(k) + Ky(k ) 
y(k) = Cx(k) + Du(k) + e(k) 


Define 


A = A- KC , B = [B-KD,K] 


v(k) = 


u{k) 
y(k ) 


(94) 


Equation (93) becomes 


x(k + 1) = Ax{k) + Bv(k) 
y(k) = Cx(k) + Du(k) + e(k) 


(95) 


Provided that p x is sufficiently large such that 

A k = {A-KC) k =0 k>p x (96) 

the input-output description with zero initial conditions can be approximated by 


y(k) - £ C{A - KCf' Ky{k - i ) = ^C(A - KC )‘~ 1 (B - KD)u(k - i ) + Du(k) + e(k) (97) 

i=l i=I 

This equation has the following linear difference structure 


y(k) + £ Ay(k -/) = £ BMk - i) + e(k) 

i=l i'=0 


(98) 


where 


A i =~C{A-KCf i K , 5o = D 

B-rrC^-ATC)'' 1 ^ , / = !, 2, p, 


(99) 


Observe that this is an extension of the deterministic linear difference equation presented in Eq. (1). 
In Eq. (98), the term e(k) tends to the residual of the Kalman filter provided that Pi is sufficiently 
large such that the condition given in Eq. (96) holds. The stochastic model structure in Eq. (98), 
which is only valid if P\ is sufficiently large, has an important interpretation in the deterministic 
case which does not require p i to be as large. This can be seen be re-examining the deterministic 
formulation in light of the stochastic consideration here. Comparing Eq. (97) to Eq. (21) reveals 
that P\—p and K = -M , provided that e(k) = 0. Thus the matrix K plays the role of — M in the 
deterministic case. The role of the matrix M in the deterministic case lias already been explained. 
One thus see that the stochastic equation considered in Eq. (97) takes a precise meaning when 
reduced to the deterministic case. Consequently, the deterministic solution for the observer 
Markov parameters given in Eq. (81) or Eq. (82) has a precise interpretation with respect to the 
Kalman filter when applied to noisy data. This analysis only serves to interpret the meaning of the 
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residual term in Eq. (98). In system identification, in order to assure that the identified residual 
indeed converges to the Kalman filter residual, the data set used must be sufficiently long. 

The Markov parameters associated with a Kalman filter are called Kalman filter Markov 
parameters. Identification of these parameters recovers not only the system but also an optimal 
Kalman filter. An important advantage of the structure in Eq. (97) is that the residual e(k) appears 
as an additive term. However, the model structure is valid in describing the input-output relation 
of the original system given in Eq. (91) in term of the Kalman filter only when p x is sufficiently 
large such that the condition given in Eq. (96) holds. Typically p x is required to be several times 
larger than the order n of the system, p x » n. 

8.2 ARMAX Model via a Kalman Filter 

In this section, a more general stochastic linear difference model is considered. This model 
has the form 


y{k) + 2 Ay(k - 0 = X B ‘ u ( k - *) + “ o ( 10 °) 

1=1 1=0 i =0 

This stochastic linear difference model, also known as an Auto-Regressive Moving Average with 
eXogeneous input model (ARMAX), is widely used in adaptive estimation, filtering, and control 
literature. An interpretation of the coefficients of this model in terms of the Kalman filter is given 
below. 

Add and subtract the term My(k) to the right hand side of the state equation in Eq. (92), 

x(k + 1) = Ax(k) + Bulk) + My(k ) - My(k) + Ke(k) 

= Ax(k) + Bulk) + M[Cx(k) + Du(k) + £(*)] - My(k) + Ke(k) (101) 

= (A + MC)x(k) + (B + MD)u(k ) - My(k) + {M + K)e(k) 

Define 

A = A + MC , B=\B + MD, -Ml , , M = M + K (102) 

L;y(*)J 

the original system in Eq. (91) can be described in the form 

x(k + 1) = A 'x(k) + B v(k) + Me(k) 

y(k) = Cx(k)+Du(k) + e(k) (103) 

The input-output description of the above system with zero initial condition is 
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( 104 ) 


y(k ) + £ C( A + MC) ,_1 Afy (* - 0 = £ C(A + MCf 1 (5 + MD)w(* - i) + Du(k) 

l'=l 1=1 

+ £c(A + Mcy-\M + K)e(k - i ) + e(k) 

i= 1 

where M and /? 2 are such that 


A* = (/l -l- AfC)* = 0 , kZp 2 (105) 

For an observable system of order n, it is possible that p 2 <n« p x . The coefficients of this 
ARMAX model in Eq. (100) given in terms of A, B, C, M, and K are 

4 = C(A + MCf 1 M , i = 1, 2, .... p 2 

Bo = D , B, = C(A + MC)‘ -1 M (106) 

Co = / , Q = C( A + MCf 1 (M + K ) 

This model has the advantage that Pi need not be as large as Pi in the previous model. However, 
this advantage is off-set by the introduction of the noise dynamics described by the C, coefficients. 
This makes the identification problem become non-linear in the parameters since both the C, 
coefficients and the residual e(k) are not known from the beginning. Note that in this description, 
the matrix M is used to limit the number of coefficients 4, B t , Q in Eq. (100) due to the term 
A + MC, including the case where A + MC is deadbeat. The matrix M does not play the role of an 
observer gain for the system given in Eq. (91), which is now given by AT. The algebraic relations 
between the coefficients A,-, fl, of the ARMAX model given in Eq. (106) and the system Markov 
parameters Y (k) = CA k ~ x B, the observer gain Markov parameters Y M (k) = CA k ~ l M are the same as 
given in Eq. (28) and Eq. (34). Furthermore, the Kalman filter gain Markov parameters 
Y K {k) = CA k ~ l K can also be computed from 


Y K (k) = Y (2 \k) + f d Y m (i)Y K (k-i) + C k 
1=1 


(107) 


where Y (2) (k) = 0 for k > p. Application of a realization procedure to the combined sequence of 
Markov parameters 


P(k) = [CA k -'B CA k ~'M CA k ~ l K\ 

= CA k -'[B M AT], * = 1, 2, ... (108) 

will realize simultaneously the matrices A, B, C, M, and K. 

As a final note, compare the ARMAX model in Eq. (100) with the model considered in Eq. 
(98). If the model in Eq. (98) is used, but Pi is not sufficiently large, then the residual in the 


28 


model will be colored. In fact, this colored residual can be modelled by the the C, coefficients in 
the ARMAX model given in Eq. (106). 


9. Concluding Remarks 


This paper presents an overview of Markov parameters with particular emphasis on their roles 
in system identification theory. In discrete-time analysis, the Markov parameter sequences have an 
important physical interpretation as they are the unit pulse response functions of a linear system. 
The close relationship between the Markov parameters of a system and its various input-output 
representations is a basis in system identification theory. In fact, a realization of a linear time- 
invariant system can be defined as a factorization of a sequence of Markov parameters into a 
specifically prescribed form. When driven by certain specialized inputs, the system sampled 
response can be conveniently interpreted as Markov parameters. Recent results in system 
identification from general inputs generated renewed interests in Markov parameters. In particular, 
the Markov parameters can be extended to include observer Markov parameters and Kalman filter 
Markov parameters. These concepts have been found to be quite useful in developing identification 
techniques for application in complex flexible structures. 
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